/*
 * Copyright (c) Jan Šlahora (${USER}) - 2011.
 * Created in IntelliJ IDEA 10.5
 */

package statistics;

import com.sun.xml.internal.bind.v2.runtime.output.SAXOutput;

import javax.security.auth.callback.TextOutputCallback;
import java.math.BigDecimal;
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import java.util.Set;

/*
*   File Pocitani 
*   in package stuff
*   was created by slaha
*   on 15.12.11 at 10:12.
*
*/
public class Pocitani {

  public double getChiKritickaHodnota(int stupne, double prav) throws Exception {
    return new Chi().chiKritickaHodnota(stupne, prav);
  }

  public double getFKritickaHodnota(int x, int y, double prav) throws Exception {
    return new FF(1-prav, x, y).getF();
  }
  public static double kombinacniCislo(int n, int r) throws Exception {
    BigInteger nF = faktorial(n);
    BigInteger rF = faktorial(r);
    BigInteger rozdilF = faktorial((n-r));

    BigDecimal citatel = new BigDecimal(nF);
    BigDecimal jmenovatel = new BigDecimal(rF.multiply(rozdilF));

    return citatel.divide(jmenovatel).doubleValue();
  }

  public static double kombinacniCislo(int n, int[] r) {
    BigInteger citatel = new BigInteger("0");
    BigInteger jmenovatel = new BigInteger("0");
    try {
      int sumaR = 0;

      BigInteger[] rF = new BigInteger[r.length];
      for (int i = 0; i < r.length; i++) {
        sumaR += r[i];
        rF[i] = faktorial(r[i]);
      }

      citatel = faktorial(n);
      jmenovatel = faktorial(n - sumaR);


      for (BigInteger bI: rF) {
        jmenovatel = jmenovatel.multiply(bI);
      }
      } catch (Exception e) {
        e.printStackTrace();
      }
    return citatel.divide(jmenovatel).doubleValue();
  }

  public static BigInteger faktorial(int x) throws Exception {
      if (x<0) {
        throw new Exception("It is not possible to calculate the factorial of a negative number - [" + x + "]");
      }

    BigInteger n = BigInteger.ONE;
    for (int i = 1; i <= x; i++) {
      n = n.multiply(BigInteger.valueOf(i));
    }
    return n;
  }

  public static int[] convertIntegers(List<Integer> integers) {
    int[] ret = new int[integers.size()];
    for (int i = 0; i < ret.length; i++) {
      ret[i] = integers.get(i).intValue();
    }
    return ret;
  }
  
  public static double chiKvadrat(Object[] cetnosti, Object[] teorCetnosti) {
    double suma = 0;

    for (int i = 0; i < cetnosti.length; i++) {
      String str = teorCetnosti[i].toString();
      double npi = new Double(str.replace(',', '.'));
      double mocnina = Math.pow((new Double(cetnosti[i].toString()) - npi), 2);
      mocnina /= npi;
      suma += mocnina;
    }
    return suma;
  }

  public static double zaokrouhleni(double cislo, int desetinnaMista) {
    if (desetinnaMista < 0) {
      return cislo;
    }

    int exp = (int) Math.pow(10, desetinnaMista);
    double d = cislo * exp;
    d = Math.round(d);
    return (d/exp);
  }

  private class  Chi {

    double BIGX = 20.0; /* max value to represent exp(x) */
    int stupne;
    double pravdepodobnost;

    private Chi () {
      super();
    }
    
    public double chiKritickaHodnota(int s, double prav) throws Exception {
      this.stupne = s;
      this.pravdepodobnost = prav;
      
      if (pravdepodobnost < 0 || pravdepodobnost > 1) {
        throw new Exception("Pravděpodobnost nemůže být menší než 0 a větší než 1!");

      }

      if (pravdepodobnost == 1) {
        return Double.MAX_VALUE;
      }

      pravdepodobnost = 1 - pravdepodobnost;

      double  CHI_EPSILON = 0.000001;   /* Accuracy of critchi approximation */
      double CHI_MAX = 99999.0;        /* Maximum chi-square value */
      double minchisq = 0.0;
      double maxchisq = CHI_MAX;
      double chisqval;

      if (pravdepodobnost <= 0.0) {
        return maxchisq;
      } else {
        if (pravdepodobnost >= 1.0) {
          return 0.0;
        }
      }

      chisqval = stupne / Math.sqrt(pravdepodobnost);    /* fair first value */
      while ((maxchisq - minchisq) > CHI_EPSILON) {
        if (pochisq(chisqval, stupne) < pravdepodobnost) {
          maxchisq = chisqval;
        } else {
          minchisq = chisqval;
        }
        chisqval = (maxchisq + minchisq) * 0.5;
      }
      return chisqval;
    }

    private double pochisq(double x, int df) {
      double a, y, s;
      double e, c, z;
      boolean even;                     /* True if df is an even number */

      final double LOG_SQRT_PI = 0.5723649429247000870717135; /* log(sqrt(pi)) */
      final double I_SQRT_PI = 0.5641895835477562869480795;   /* 1 / sqrt(pi) */

      if (x <= 0.0 || df < 1) {
        return 1.0;
      }

      a = 0.5 * x;
      even = (df % 2 == 0)? true : false;
      y = 0;
      if (df > 1) {
        y = ex(-a);
      }
      s = (even ? y : (2.0 * poz(-Math.sqrt(x))));
      if (df > 2) {
        x = 0.5 * (df - 1.0);
        z = (even ? 1.0 : 0.5);
        if (a > BIGX) {
          e = (even ? 0.0 : LOG_SQRT_PI);
          c = Math.log(a);
          while (z <= x) {
            e = Math.log(z) + e;
            s += ex(c * z - a - e);
            z += 1.0;
          }
          return s;
        } else {
          e = (even ? 1.0 : (I_SQRT_PI / Math.sqrt(a)));
          c = 0.0;
          while (z <= x) {
            e = e * (a / z);
            c = c + e;
            z += 1.0;
          }
          return c * y + s;
        }
      } else {
        return s;
      }
    }

    private double ex(double x) {
      return (x < -BIGX) ? 0.0 : Math.exp(x);
    }

    private double poz(double z) {
      double y, x, w;
      double Z_MAX = 6.0;              /* Maximum meaningful z value */

      if (z == 0.0) {
        x = 0.0;
      } else {
        y = 0.5 * Math.abs(z);
        if (y >= (Z_MAX * 0.5)) {
          x = 1.0;
        } else if (y < 1.0) {
          w = y * y;
          x = ((((((((0.000124818987 * w
              - 0.001075204047) * w + 0.005198775019) * w
              - 0.019198292004) * w + 0.059054035642) * w
              - 0.151968751364) * w + 0.319152932694) * w
              - 0.531923007300) * w + 0.797884560593) * y * 2.0;
        } else {
          y -= 2.0;
          x = (((((((((((((-0.000045255659 * y
              + 0.000152529290) * y - 0.000019538132) * y
              - 0.000676904986) * y + 0.001390604284) * y
              - 0.000794620820) * y - 0.002034254874) * y
              + 0.006549791214) * y - 0.010557625006) * y
              + 0.011630447319) * y - 0.009279453341) * y
              + 0.005353579108) * y - 0.002141268741) * y
              + 0.000535310849) * y + 0.999936657524;
        }
      }
      return z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5);
    }
  }

  private class FF {
    //http://faculty.vassar.edu/lowry/tabs.html

    double F;

    public double getF() {
      return F;
    }

    FF(double p,int n1, int n2) {
      double v=0.5;
      double dv=0.5;
      double f=0;

      while(dv>1e-6) {
        f=1/v-1;
        dv=dv/2;
        if(fishF(f, n1, n2)>p) {
          v=v-dv;
        } else {
          v=v+dv;
        }
      }
      this.F = f;
    }

    private double fishF(double f,int n1, int n2) {
      double x=n2/(n1*f+n2);

      if((n1%2)==0) {
        return StatCom(1-x,n2,n1+n2-4,n2-2)*Math.pow(x,n2/2);
      }

      if((n2%2)==0) {
        return 1-StatCom(x,n1,n1+n2-4,n1-2)*Math.pow(1-x,n1/2);
      }

      double PiD2=Math.PI/2;
      double th=Math.atan(Math.sqrt(n1*f/n2));
      double a=th/PiD2;
      double sth=Math.sin(th);
      double cth=Math.cos(th);

      if(n2>1) {
        a=a+sth*cth*StatCom(cth*cth,2,n2-3,-1)/PiD2;
      }

      if(n1==1) {
        return 1-a;
      }

      double c=4*StatCom(sth*sth,n2+1,n1+n2-4,n2-2)*sth*Math.pow(cth,n2)/Math.PI;

      if(n2==1) {
        return 1-a+c/2;
      }

      double k=2;

      while (k<=(n2-1)/2) {
        c=c*k/(k-.5); k=k+1;
      }

      return 1-a+c;
    }

    private double StatCom(double q, int i,int j, int b) {
      double zz=1;
      double z=zz;
      int k=i;

      while(k<=j) {
        zz=zz*q*k/(k-b);
        z=z+zz; k=k+2;
      }
      return z;
    }
    
  }

  public static void main(String[] args) {
    Pocitani p = new Pocitani();


 try {

    System.out.println(p.getChiKritickaHodnota(4, 0.95));
   System.out.println(p.getFKritickaHodnota(1, 3, .05));
   System.out.println(p.getFKritickaHodnota(4, 30, .05));
    } catch (Exception e) {
      e.printStackTrace();  //To change body of catch statement use File | Settings | File Templates.
    }
  }
}
